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We derive the exchange-correlation potential corresponding to the nonlocal van der Waals density 
functional [M. Dion, H. Rydberg, E. Schroder, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. 
Lett. 92, 246401 (2004)]. We use this potential for a self-consistent calculation of the ground state 
properties of a number of van der Waals complexes as well as crystalline silicon. For the latter, where 
little or no van der Waals interaction is expected, we find that the results are mostly determined 
by semilocal exchange and correlation as in standard generalized gradient approximations (GGA), 
with the fully nonlocal term giving little effect. On the other hand, our results for the van der Waals 
complexes show that the self-consistency has little effect at equilibrium separations. This finding 
validates previous calculations with the same functional that treated the fully nonlocal term as a post 
GGA perturbation. A comparison of our results with wave-function calculations demonstrates the 
usefulness of our approach. The exchange-correlation potential also allows us to calculate Hellmann- 
Feynman forces, hence providing the means for efficient geometry relaxations as well as unleashing 
the potential use of other standard techniques that depend on the self-consistent charge distribution. 
The nature of the van der Waals bond is discussed in terms of the self-consistent bonding charge. 

PACS numbers: 31.15.Ew, 71. 15. Mb, 61.50.Lt 
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I. INTRODUCTION 

Density functional theory (DFT) has acquired high 
respect for its simplicity and accuracy within first- 
principles calculations. In particular, generalized gradi- 
ent approximations^^ have had much success in describ- 
ing isolated molecules^ as well as dense bulk matter £ For 
molecules, this performance has been further improved 
by the use of various empirical and hybrid methods. On 
the other hand, for van der Waals (vdW) complexes and 
sparse matter, including many layered structures, poly- 
mer crystals, and organic molecular crystals, these meth- 
ods either give sporadic results or fail completely. The 
reason for this failure is that the dominant part of the 
stabilization energy in many cases comes from the dis- 
persion energy, which is not correctly accounted for in 
standard DFT&2£&±£ 

To remedy the situation, a new approach using a van 
der Waals density functional (vdW-DF) with a nonlo- 
cal correlation energy has been developed* 1 ^ This for- 
malism includes van der Waals forces in a seamless fash- 
ion and has been applied with quite good results to 
bulk layered systems^ dimerS ) 10 ' 11 : 12 ' 13 : 14 ' 15 : 16 molecules 
physisorbed on infinite surfaces^ and bulk crystals of 
polyethylene* 1 ^ Lacking the exchange-correlation poten- 
tial corresponding to the vdW-DF, these calculations 
evaluated the vdW-DF as a post-processing perturba- 
tion, using a charge density obtained from self-consistent 
calculations with the PBE or revPBE functionals. 2,3 

The goal of this paper is to derive the exchange- 
correlation potential for the vdW-DF, thereby enabling 
fully self-consistent calculations using vdW-DF. This 



exchange-correlation potential will give more credibil- 
ity to the post-process vdW-DF calculations mentioned 
above. Furthermore, it encourages us to apply the vdW- 
DF to even more exciting cases in the future. In addi- 
tion, the knowledge and implementation of the exchange- 
correlation potential provides the missing underlying 
framework for further developments, such as the eval- 
uation of forces and the stress tensor. These tools are 
of tremendous interest for the structural optimization of 
many soft materials. Also many other standard tools 
for analyzing materials properties that are usually built 
into ab-initio codes become accessible through the self- 
consistent vdW-DF potential. With these developments, 
vdW-DF should become a standard tool for systems 
where vdW interactions are important. This includes not 
only extended systems, but also finite systems that are 
too large for the standard quantum chemical methods to 
be effective. 

For such large finite systems, the advantage of the 
vdW-DF approach becomes obvious when the scaling of 
the computational effort with respect to system size is 
considered. Particular interest in van der Waals com- 
plexes has emerged in the context of large systems, 
where questions regarding DNA base-pair bonding, pro- 
tein structure and folding, and organic molecule crystal- 
lization are experiencing a surge of interest. Tradition- 
ally, second-order M0ller-Plesset perturbation theory^ 
(MP2) and coupled-cluster calculations^ with singles, 
doubles, and perturbative triple excitations (CCSD(T)) 
would be used to study these systems as they are re- 
garded as accurate and reliable. Unfortunately, in many 
of these cases their application is limited to relatively 
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small systems, since they are too computationally de- 
manding. Furthermore, the results are basis set depen- 
dent and elaborate techniques are necessary to estimate 
results at the complete basis set limit. This is where the 
strength of the vdW-DF approach lies: Since it scales 
with system size N just like a regular DFT calculation — 
usually 0(N 3 ), complexes currently out of the reach for 
MP2 and CCSD(T) calculations are easily treatable. 

For extended systems that require the consideration of 
long-range van der Waals interactions, the advantage of 
vdW-DF is even more obvious. To our knowledge there 
is simply no other first-principles method available. The 
earlier post-process application of the method has al- 
ready given reasonable results for a polymer crystal* 
and layered intercalates^ and applications to molecular 
crystals are under way. The results of the this work sug- 
gest that the predictions will not be greatly affected by 
the full self-consistency that now becomes possible. 

The vdW-DF treats the nonlocal vdW interaction by 
an approximate non-empirical method that is not only 
correct in principle, but useful in practice. It differs from 
other methods^ 3 -^^^^^^ that treat the vdW in- 
teraction as occurring directly between the nuclei via em- 
pirically determined potentials. As discussed later, the 
correct source for the forces on the nuclei is electrostatic, 
and results from a change in the charge densities of the 
fragments as the vdW bond is formed. This effect is in- 
cluded correctly in vdW-DF. 

We have organized this paper in the following way. In 
Sec. [II] we analytically derive the exchange-correlation po- 
tential. In Sec. lIIII the potential is then applied to several, 
very different physical systems. We compare our results 
with other calculations such as MP2 and CCSD(T) and 
find good agreement. Our results show that calculating 
vdW-DF binding energies in a post-process procedure, 
rather than self-consistent, is reasonable, thus justifying 
the approximation used in earlier calculations. Here, we 
also deal with the question of how the vdW-DF performs 
for solids and covalently bonded systems such as a CO2 
molecule or crystalline Si. In Sec. IIVI we use the vdW- 
DF self-consistent charge density to discuss the nature of 
the van der Waals bond. We conclude in Sec. [Vj Some 
details of the derivation of the potential are collected in 
Appendices [£] |Bj and [C] 



Waals interaction enters through a fully nonlocal correc- 
tion [n] that concerns the correlation only. In the 
following we are interested in the exchange-correlation 
potential that originates from this nonlocal contribu- 
tion. However, we have to keep in mind that the total 
exchange-correlation potential corresponding to Eq. ([T]) 
also includes the parts stemming from £'JJ ovPBE [n] and 
E^ BA [n] — these expressions are well-known and can be 
found elsewherej 30 ' 31 

In order to calculate the potential of interest, we have 
to take the functional derivative of the energy with re- 
spect to the density, i.e. 
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where 4>(r, r') is a function depending on r — r' and the 
electronic densities n and the magnitude of their gradi- 
ents at the points r and r'. Details about the kernel 
(j)(r, r') and its evaluation can be found in Appendix [Al 
It follows that 



vf(r) = = / d 3 rd 



5n(r) 



(r,r>(r') 



-I — / d rd r n(r)<£(r, r ) 



Sn(r') 



5n(r) 



(4) 



on(r) 



The first two lines can be simplified if we use 
6n(r) 
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It turns out that the kernel </>(r,r') depends on r and r' 
only through two functions d and d' 
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where we have introduced R rr / = r— r' and R rr i — |R rr ' |. 
The definition of qo(r) can be found in Eqs. (11) and (12) 
of Ref. 10] which reads 



II. DERIVATION OF THE 
EXCHANGE-CORRELATION POTENTIAL 

The density functional in question is defined in 
Ref. [lOj. It consists of several different contributions, 
i.e. 
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The first two parts are simply revPBE exchange and LDA 
correlation. We note that these parts represent LDA plus 
gradient corrections and are hence semilocal. The van der 
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where we have used the standard expressions for the 
Fermi wave vector and the reduced gradient 
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and Z a t = —0.8491. A complete discussion of the gradi- 
ent term in Eq. ([7]) is given in Appendix IB1 
We may now write 
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The symmetry in the last equality is apparent by exploit- 
ing the definition of cf)(r, r') in Appendix lAl In order to 
calculate the derivative in the third line of Eq. (j4]) we can 
use this symmetry and find 
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If we keep in mind that both terms appear under a dou- 
ble integral, we can change the integration variables and 
find that these two terms are equivalent. Then, sim- 
ply inserting a factor of two and abbreviating d ^g^ d - = 
<pd(d,d r ) = 0d(r,r') yields 
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The missing derivative in the equation above follows 
from the definition of d(r,r') in Eq. ([6]), i.e. 
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If we denote the derivative of s]^ A with respect to the 
density as £x,P A an d simplify further, it follows 
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Equation (fT"3f together with Eq. (|T2|) can now be inserted 
in the expression for the potential, Eq. (fTTj) . and it follows 
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The third line in this equation includes the gradient of 
the 5-function. Partial integration lets us rewrite the 
corresponding integral as 
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where the sufficiently rapid fall- off of <f> d (see Eq. l2"0")) guar- 
anties the vanishing of the integrated part. All contribu- 
tions now have the common form J d 3 r' . . . n(r') so that, 



after dropping the tilde (r — > r), we may simply write: 
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The term including the gradient V •[...] in the equa- 
tion above can be further simplified, but the details are 
deferred to Appendix O If we use Eq. © to replace R rr > 
in favor of d(r,r') or d'(r, r') and collect similar terms, 
we find as the final result 
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where the functions an (r, r') and $i(r,r') are given by: 
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where $4 is the kernel function defined in Ref. [lOj. Here, 
we have abbreviated higher order derivatives of <p(d,d') 
as previously, 

9 2 0(d,d') d 2 0(d,d') 

= ^ dd (d, d ) and = (?W(d, d ) . 
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For completeness, details about the evaluation of 
4> d {d,d'), <j)dd(d,d'), and 4> d d'(d,d') are provided in Ap- 
pendix lAl 

The three extra internal kernel functions <&i (d, d') 
through $3(d, d') in Eqs. (|18aH18c]l necessary for f" 1 are 
the analogues of the single function <p(d, d') (= $4(d, d')) 
used for However, unlike the latter, they are not 

symmetric under the exchange of d and d'. Consequently, 
when expressed in terms of the sum and difference vari- 
ables D and S defined by-i£ 



d = D(l + S) 
d! = D(l-S) 
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FIG. 1: The kernels $1, $2, and $3 given in Eqs. (fl8|) are 
plotted as a function of D for values of 5 = ±0.9, ±0.5, and 
0.0 [see Eqs. 1(15)1 ]. 



their values depend on the sign of <5. Plots of these three 
new kernels are given in Fig. [TJ Comparison with Fig. 1 
of Ref. [HI (see erratum), shows that there are no unex- 
pected features. The details of the evaluation are given 
in Appendix [AJ The asymptotic form of the kernels for 
w" 1 for large R = |r — r'| may be determined by using the 
relation^ 
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The use of Eq. ([20]) in Eqs. CEH]) plus the fact that D cx R 1 
and S oc R° gives the asymptotic dependence R~ 6 to $1, 
$2, and $4, while $3 falls off faster and thus does not 
contribute to the asymptote. 

The asymptotic form above follows from the form of 
the model response function that determines the kernels. 
It gives the correct asymptotic dependence for all finite 
sized fragments, and for all infinite fragments that are 
not metallic. A couple of decades ago, Barash^S, showed 
that infinite metallic fragments of reduced dimensional- 
ity could show a more slowly decaying asymptote than 
would be given by the appropriate integrals over (|20[) . 
This effect occurs because the dipole response functions 
of infinite metals may not be sufficiently convergent at 
small frequency and wave vector. A number of examples 
of special cases have been worked out more recently^ We 
can say with certainty that, with the response function 
presently used, our functional will fail for infinite metal- 
lic (or semimetallic) fragments of reduced dimensional- 
ity at very large distances (for example, parallel metallic 
sheets). The distance where the crossover occurs will de- 
pend on the details of the material. We are unaware of 
any such material where this distance is known. 



The arrangement of the results for w" 1 in Eq. (|T6|) ac- 
cording to Eqs. (fTT)) and (fT%)) was chosen to facilitate 
rapid numerical evaluation. For a grid of N points, the 
evaluation of u" 1 on each is 0(N 2 ) according to Eq. (|16|) . 
However, the evaluation of the <&'s in Eqs. (|18[) needs only 
to be done once to create an interpolation table. The 
conversion between r (r') and d (d') [Eqs. ©] of course 
needs to be done 0(N 2 ) times, but the lengthy part of 
that would be the calculation of qa, which instead can be 
calculated in advance on each grid point, i.e. it is O(N). 
The calculation of ati and a 2 [Eqs. JT7J)] is O(N). So, 
the only parts of the calculation that are 0(N 2 ) are the 
calculation of the components of R for a.3 [Eq. (|17c[) ]. 
the calculation of |R| for r to d conversion [Eqs. ©], two 
table interpolations, and a few arithmetic operations. 



III. RESULTS 

The exchange-correlation potential derived in the pre- 
vious section will now be applied to various systems. Cal- 
culations labeled as self-consistent (SC) include the po- 
tential as part of the Kohn-Sham potential, so that 
the density used to evaluate the nonlocal correlation en- 
ergy, Eq. is fully self-consistent. Results in which 
the nonlocal vdW correlation energy ([2]) is calculated us- 
ing the density from a self-consistent PBE calculation 
will be referred to as non-SC. It should be noted that 
all previous vdW-DF calculations^^^±^I^>i&±^ 
were preformed with such a non-SC post-process proce- 
dure. 

It is apparent that the post-process evaluation of the 
vdW-DF is an approximation since the charge density is 
not allowed to evolve under the full vdW-DF functional. 
Although there are many arguments that the effect of this 
approximation is small, the ultimate test is obviously a 
comparison of our new SC results with the older non-SC 
results. In the following we shall make this comparison. 

In addition to the SC and non-SC results, we will also 
show results marked as "no Ef". These results were 



obtained by only considering E" 



rcvPBE 



E. 



LDA 



Eq. © 



and neglecting the nonlocal part E~. This will allow us 
to prove that the correct behavior of the vdW-DF for 
van der Waals systems indeed is encoded in the nonlocal 
correlation E^ 1 . The calculations below were made with 
the aid of the plane-wave code ABINIT 34 and the real- 
space code PARSEC^ 



A. Ar and Kr dimers 

We start out by applying the self-consistent vdW-DF 
to some rare gas dimers, which, due to their closed va- 
lence shells, owe much of their binding to van der Waals 
interactions. Results for the interaction energy of an Ar 
and a Kr dimer as a function of separation are depicted in 
Figs. and [3J It can be seen that the differences between 
SC and non-SC calculations are indeed very small, and 
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FIG. 2: (top) Interaction energy of the Ar dimer as a func- 
tion of separation. Plotted are self-consistent and non self- 
consistent results. In addition, we show results where has 
been neglected. Experimental values are taken from Ref. [3(|. 
(bottom) Forces calculated as the derivative of the energy 
(—dE/dr) and the Hellmann-Feynman forces. 



on the scale of this plot they are negligible. In general, 
these differences are larger for smaller separations, where 
the perturbation from the nonlocal term is larger. For 
Ar and Kr we find the binding distances to be 3.9 A and 
4.2 A, which is approximately 5% larger than their corre- 
sponding experimental values^ On the other hand, our 
interaction energies are stronger than experiment. This 
behavior of the binding distance is consistent with all 
previous calculations^ ' 10 ' 11 ' 12 ' 13 ' 14 ' 15 ' 16 ! 17 ! 18 and has been 
attributed to the form of the exchange usedj 12 ' 13 In ad- 
dition, it can be seen that without the nonlocal contri- 
bution (no E™ 1 ) the binding is much too weak and the 
separation is too large. In some cases, as we shall see later 
in the text, without this contribution binding does not 
occur at all. It thus becomes obvious that the addition 
of the nonlocal contribution i?" 1 dramatically improves 
the description of these van der Waals systems. 

As a result of our DFT calculations, we not only obtain 
the energies discussed above, but also the correspond- 
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FIG. 3: Same as Fig. [2] except here describing the Kr dimer. 



ing Hellmann-Feynman forces which act on each nuclei. 
These are simply the electrostatic forces, which require 
the correct self-consistent charge distribution for their de- 
termination. Until now, it was impossible to obtain them 
within vdW-DF. Results for these forces are plotted in 
the bottom panels of Figs. [2] and [3] In order to verify the 
correctness of these forces, we have plotted the negative 
derivative of the energy from above (—dE/dr). It can be 
seen that the zero point of the forces perfectly aligns with 
the corresponding energy minimum. As for the energy, 
we have again included results for "no . It is now 
apparent that semilocal exchange-correlation alone can- 
not describe van der Waals systems correctly and most 
of the binding force originates in the nonlocal part of the 
functional. 



B. C0 2 dimer 

Next, we apply the vdW-DF to a slightly more com- 
plex molecule, i.e. the CO2 dimer. However, before we 
discuss the details about the dimer, it is interesting to 
first explore the behavior of the vdW-DF for a single 
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FIG. 4: Interaction energy of a CO2 dimer as a function of its 
separation at a fixed slip distance of 2.05 A. The inset shows 
the CO2 dimer geometry. 



CO2 molecule. The CO2 molecule is a covalently bonded 
system, with only minor contributions from vdW interac- 
tions. As such, it is interesting to investigate how vdW- 
DF treats such a system. The important question here is 
whether or not the functional will continue to give good 
results for cases in which van der Waals interactions are 
negligible or non-existent. 

For CO2 we can partly answer this question by calcu- 
lating the optimal C-0 bond length. Our results show 
that the bond length that vdW-DF predicts is almost 
identical to the bond length obtained from a standard 
PBE calculation. It turns out that the nonlocal part of 
the functional only makes a negligible contribution 

to the total energy. In this regard, the total exchange- 
correlation energy computed with the vdW-DF is just 
the sum of the £™ vPBE [n] and E% DA [n], i.e. revPBE ex- 
change and LDA correlation. It remains to be determined 
whether this form of exchange-correlation energy is ap- 
propriate for the system of interest, but the main point is 
that in the regime of weak van der Waals forces, the vdW- 
DF will not change or alter well established DFT results. 
We will readdress the question of covalently bonded sys- 
tems and vdW-DF in the context of a solid in Sec. IIII El 

Now we turn to the bonding of the CO2 dimer. We 
have performed calculations for the interaction energy 
of this dimer as a function of the dimer separation and 
slip distance. We find that the preferred parallel slip dis- 
tance is 2.05 A, and the corresponding interaction energy 
as a function of dimer separation is depicted in Fig. 2] 
Similar to Ar and Kr, the results indicate that the differ- 
ences between non-SC and SC vdW-DF calculations are 
negligible. Furthermore, the C-C separation distance of 
3.772 A is in good agreement with the experimentally 
determined separation distance of 3.602 A.— The vdW- 
DF interaction energy of 1.55 kcal/mol is comparable 
to the MP2 value of 1.36 kcal/mol (extrapolated to the 
complete basis set limit )2S as well as the value of 1.60 
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FIG. 5: Interaction energy of water on benzene. The little 
inset shows the geometry. 

kcal/mol determined by Becke and coworkers using their 
density-functional model for the dispersion interaction^ 
In addition, Fig. [4] shows that the equilibrium separation 
distance for the CO2 dimer is slightly larger than the 
MP 2 results^ 

As before, by looking at the "no E 1 " 1 " curve we can see 
that most of the binding originates in the nonlocal part 
of the functional. 



C. Water on benzene 

We now move to the physically and chemically more 
interesting question of how a typical solvent like water 
binds to hydrocarbons. Here, we consider the simple test 
case of the interaction energy between water and ben- 
zene. We limit ourselves to a particular configuration, 
in which the water molecule is positioned with the hy- 
drogen atoms at equivalent heights pointing towards the 
benzene plane. The oxygen atom is positioned directly 
above the center of the benzene molecule, as depicted in 
the little inset in Fig. [51 

First, we relaxed the bond length of the water molecule 
and benzene ring individually. The optimal H-0 bond 
length for the water molecule and the H-O-H angle were 
found to be 0.973 A and 106°, respectively, in quite good 
agreement with the experimental values of 0.958 A and 
104.5 . 40 The C-C and C-H bond lengths for the benzene 
molecule were 1.396 A and 1.100 A, where experiments 
find 1.397 A and 1.084 A, respectively^ Our results for 
the interaction energy versus the vertical separation dis- 
tance between the oxygen atom and the center of the 
benzene ring is plotted in Fig. [5] Again, we see that the 
difference between SC and non-SC calculations is very 
small. Consistent with our previous results, the differ- 
ence is more evident at shorter separation distances. For 
comparison, we plotted the results from a recent MP2 
studyi^! We find fair agreement, with MP2 calculations 
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FIG. 6: Interaction energy of a cytosine dimer as a function 
of separation. 



giving a shorter equilibrium distance and lower interac- 
tion energies than vdW-DF. The difference for both is 
less than 5%. 

Once again, we find that without E£ , the complex 
binds at a much too large distance with much weaker 
binding energies — supporting the fact that the nonlocal 
correlation energy is crucial for the correct description of 
van der Waals binding. 



D. Cytosine dimer 

The applicability of vdW-DF to monosubstituted ben- 
zene dimers has been shown in Rcf. [13]. We will now 
extend this work by considering cytosine, a simple nu- 
cleic acid base. 

We study the interaction energy of a cytosine dimer as 
a function of separation (Fig. [6]). The geometry is such 
that the two molecules are placed on top of each other 
with one of them rotated by 180° around an axis that 
passes through the center of both rings. Consistent with 
other findings reported in this paper, there are no sig- 
nificant differences between the SC and non-SC results. 
Similarly as in the cases before, the "no -E" 1 " curve re- 
veals that the majority of the binding energy comes from 
the nonlocal contribution of the vdW-DF. As usual, while 
the comparison to MP2 calculations is quite good, we find 
a slightly larger binding distance. 



E. Crystalline Silicon 

We now choose to apply the new functional to an 
extended system and return to the question raised in 
Sec. IIIIB1 We purposefully selected an extended system 
in which the bonding is mostly covalent and where van 
der Waals interactions should be unimportant. Again, 




4.5 5.0 5.5 ^ 6.0 6.5 

Lattice parameter [A] 



FIG. 7: Total energy (per two-atom unit cell) of Si as a func- 
tion of the lattice constant. The inset shows the convergence 
of with respect to the cutoff radius i? m ax. See text for 
further details. 



we seek to understand how the vdW-DF treats such sys- 
tems. 

In Fig. [71 we plot the total energy (per two-atom unit 
cell) of crystalline silicon as a function of the lattice pa- 
rameter. In addition to the "usual" curves for SC and 
non-SC results, which show no differences, we now in- 
clude results from calculations performed with standard 
LDA and PBE functionals. The main point of Fig. [7J 
is that all calculations result in very similar lattice con- 
stants: non-SC (5.48 A), SC (5.48 A), no Ef (5.48 A), 
PBE (5.46 A), and LDA (5.37 A). In other words, the 
nonlocal contribution in systems like silicon is indeed neg- 
ligible and the functional performs as expected for the 
local and semilocal contributions, i.e. revPBE exchange 
and LDA correlation. This is even more evident in the 
fact that the lines for SC and "no E"" 1 " almost coincide, 
resulting in exactly the same lattice constant. 

At this point, we would also like to mention another 
important aspect of periodic systems. Since crystalline 
Si is an extended system, interactions between different 
unit cells have to be included when the energy in Eq. (jTJ) 
and the potential in Eq. (fTB]) are evaluated. From Fig. Q] 
it is apparent that all kernel functions decrease in ab- 
solute magnitude with increasing D, and in turn, also 
with increasing separation |r — r'| of the two points r and 
r'. Thus, it is natural to introduce a maximum distance 
-Rmax above which the contribution of a given pair r and 
r' is no longer included. For the case of Si, the con- 
vergence of with respect to this cutoff radius i? max 
is depicted in the inset in Fig. [7J It can be seen that 
the convergence is quicker than for materials containing 
carbonj 17 ! 18 For the purpose of this study we have chosen 

Rmax = 15 A. 

The findings for crystalline silicon parallel the earlier 
findings for the CO2 molecule in Sec. IIIIBI When the 
van der Waals interactions in an extended system be- 
come negligible, the contribution of the nonlocal part E^r 
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to the total energy becomes negligible. Again, vdW-DF 
does not affect well established DFT results. With this, 
we conclude our quantitative results and move on to a 
qualitative description of the van der Waals bond. 



IV. NATURE OF THE VAN DER WAALS BOND 

The conventional picture of the van der Waals inter- 
action envisions in its simplest form two fragments with 
distinct charge distributions, whose motions are corre- 
lated so as to reduce the interfragment electron-electron 
repulsion and hence produce a net attractive force be- 
tween the fragments. For distant fragments the dipolar 
interactions caused by these correlations are dominant, 
and one obtains the familiar attractive r~ 6 interaction 
(at very large distances, typically 100 A or more, this in- 
teraction is further weakened by relativistic retardation 
effects). The vdW-DF method embraces this picture, 
which indeed was central to its derivation, and also in- 
cludes the contribution of more complex correlated mo- 
tions that occur as the fragments get closer and even 
merge. 

However, the van der Waals bond must also have a 
completely different feature as part of its nature. This 
feature arises because the nuclei are classical particles 
influenced only by Coulomb forces and immune to the 
fluctuations in such forces due to the more rapid elec- 
tronic motions. Because of the stationary property of 
the energy with respect to variations in the wave func- 
tion, these forces are easily calculated — a result often at- 
tributed to Hellmann^A and Feynman,— but whose foun- 
dation is much olderj 43 i 44 The consequence is that the 
static or ground state electronic charge distribution must 
deform itself in such a way as to produce the required 
forces on the nuclei by classical Coulomb interactions 
alone. This concept is familiar for covalent bonds, but 
how is it implemented for van der Waals bonds? 

We can now answer this question by calculating how 
the static charge density changes when the nonlocal con- 
tribution E% 1 is included in the functional. To this end, 
we investigate the density of the Ar dimer in a plane 
that includes both atoms. We calculate the induced 
electron density that occurs when the atoms bond, i.e. 
?i ind = ribond - "atom, where n bon d is the density of the 
dimer and n at0 m is the density of the isolated atoms. 
We can now study the difference in the induced electron 
density when the full functional is used nf" d , compared 
to the induced density that results without the nonlo- 
cal part n"^ d nl . Results for n[jj d — n^ d nl for a separation 
of 4.02 A are depicted in Fig. [8] The plot shows the 
change in the induced electron density when the nonlo- 
cal part of the functional is "turned on." As expected, 
it can be seen that the electron density from around the 
nuclei moves between the atoms, which in turn explains 
the stronger binding in terms of implied changes in the 
electrostatic forces on the nuclei arising from this charge 
redistribution. Although the change in density is very 
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FIG. 8: The bonding charge of a van der Waals com- 
plex. Shown is the difference in induced electron density 
n tnd ~ n ind nl f° r the Ar dimer with a separation of 4.02 A. 
The scale is in A and the black dots mark the position of the 
nuclei. The zero level is marked by the black contour. Red ar- 
eas represent areas of electron density gain when the nonlocal 
part is included; conversely, blue areas indicate loss of elec- 
tron density. Increments between contour lines are 5x10 
electrons/ A 3 . 

small (~ 10 -4 electrons/ A 3 ), the resulting binding en- 
ergies and forces change considerably, as seen in Figs. [2] 
and [3] This is the van der Waals analogue of the bond- 
ing charge familiar from the theory of the covalent bond, 
although it is much weaker and of different shape. 



V. CONCLUSIONS 

The long-range part of the Kohn-Sham exchange corre- 
lation potential derived here allows the vdW-DF theory 
of Dion et al^ to be applied in a fully self-consistent 
manner. In previous work the vdW-DF functional had 
been applied as a post-GGA perturbation. It was ar- 
gued that this was a reasonable approximation, because 
one does not expect that the rather weak and diffuse 
vdW interaction should substantially change the elec- 
tronic charge distribution. The present work, through 
the application of the fully self-consistent theory to a 
number of van der Waals complexes, shows that this ar- 
gument was correct. In all the complexes studied, the 
predictions of the fully self-consistent theory were nearly 
indistinguishable from those obtained via post-GGA per- 
turbation. This was true near the equilibrium separa- 
tion between the fragments. However, inaccuracies in 
the perturbation method began to become noticeable 
when the fragments were pushed together more closely 
and the perturbation becomes larger. The above con- 
clusions additionally validate the many previously per- 
formed post-GGA perturbative calculations, and suggest 
that the more efficient post-GGA perturbative method 
will remain an effective tool. 

Nevertheless, the availability of a self-consistent the- 
ory is important, because it unleashes the possibility of 
applying a number of important and standard DFT tech- 
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niques whose availability depends on this self-consistency. 
The most obvious of these is the ability to calculate inter- 
nuclear forces via electrostatics. Indeed we show that the 
Hellmann-Feynman forces are given accurately in vdW- 
DF, a feature that will allow for efficient nuclear relax- 
ation methods to be employed to determine optimal ge- 
ometric structures. 

Finally, our applications to the geometries of isolated 
molecules as well as bulk silicon suggest that vdW-DF is 
really a general density functional. It gives good results 
for systems or parts thereof where van der Waals forces 
are substantial, but does not spoil the results of ordinary 
LDA-GGA type functionals when van der Waals forces 
are unimportant. 

These studies have shown the strengths of the ex- 
isting functional. However, areas for future improve- 
ments are also important. One of these is the con- 
sistent small overestimation of the separation between 
vdW bound fragments. This occurs in all the vdW com- 
plexes considered here, as well as those considered in 
previous work i 12 i 13 i 17 ' 18 The reason for this persistent 
shift has been attributed to inadequacies of the exchange 
functional ) 12 ! 13 ' 18 but a systematic study has yet to be 
done. It seems likely that the source of this discrepancy 
can be unambiguously identified and corrected in future 
work. 
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APPENDIX A: DETAILS OF THE EVALUATION 
OF THE KERNEL 



and the function T is given by 



T(w,x,y,z) = 



1 



1 



w + x y + z 
1 



(A3) 



1 



(w + y)(x + z) (w + z)(y + x) 
Furthermore, we have used the definitions 



and 



v(u) = u 2 /2h(u/d) , 
u'(u) = u 2 /2h(u/d') 



h(t) = 1 - cxp (-47rf 2 /9) 



(A4a) 
(A4b) 



(A5) 



Since cf> according to Eq. (|A1[) depends on d only via v 
and on d! only via v', it is straightforward to perform the 
required derivatives indicated in Eqs. (|18|) analytically, 
leaving only the double integral over a and b in Eq. (|A1[) 
to be done numerically. The required partial derivatives 
of T (u(a), v(b),v'(a), v'{b)) in Eq. are T d , T dd , and 
T d d' , where v and v' depend implicitly on d and d' re- 
spectively according to Eqs. (|A4[) . Thus, T d , Tdd, and 
T d d' may be expressed in terms of partial derivatives of 
T(w,x,y, z) in Eq. (|A3[) with respect to w, x, y, z, and 
the partial derivatives of v and v' in Eqs. (IA4p with re- 
spect to d and d! . This procedure gives 



Tdd 



T w v d {a) + T x v d {b) , (A6a) 
T ww u 2 (a) + T w v d d{a) + T xx v 2 d {b) 
+ T x u d d{b) + 2T wx v d (a)v d (b) , (A6b) 
T wy v' dl {a)vd{a) + T wz i> d ,{b)v d {a) 
+ T xy v' d ,{a)v d [b)+T xz v' d ,{b)v d [b) . (A6c) 



The algebraic evaluation of the derivatives in Eqs. (|A6[) 
is most simply done with a symbolic algebra program. 
Then, altering Eq. (|A1[) by the replacement of T with 
dT d using Eq. (|A6a[) gives $i (Eq. (|18ajl ). An analogous 
replacement using d T dd gives $ 2 (Eq. CHE}), while for 
$3 (Eq. (|18cp ) the replacement should be made with T d + 
dTdd + d'Tdd'- 



The expression for the kernel <p can be written as 10 

0(d,d') = / a 2 da b 2 dbW(a,b) 
Jo Jo 

x T (v(a), v(b), v'(a), v'(b)) , 



where the function W(a,b) is defined as 
W(a, b) 



a 3 b 3 

(3 — a 2 )b cos b sin a + (3 — b 2 )a cos a sin b 
+ (a 2 + b 2 — 3) sin a sin b — 3ab cos a cos b 



(Al) 



(A2) 



APPENDIX B: GRADIENT CORRECTION 

As discussed in Ref . [l(| , the long range part of the cor- 
relation energy is calculated assuming a simple single 
pole model for the necessary response function, with pa- 
rameters determined by sum rules plus the requirement 
that the correct energy density should be obtained when 
it is applied locally, i.e., according to LDA with the ap- 
propriate gradient correction. In particular, the gradient 
terms that represent an (inappropriate) attempt to ex- 
pand the van der Waals interaction in a gradient series 
should be omitted. 
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The leading gradient contribution to the exchange- 
correlation energy density may be written as 



£grad 



,e 2 kF ( Vn 



12tt \2k F i 



(Bl) 



This defines the quantity Z according to the usage of Lan- 
greth and Perde w 45 ' 46 and Langreth and Vosko.— Rasolt 
and Geldar t 48 ' 49 define a Z which is larger by an addi- 
tive constant of 4/3, which would cause a corresponding 
change in Eq. (|Bip . The sometimes confusing aspects 
of this notation have been clarified more recently^ The 
value of Z was first inferred by combining the calculations 
of Ma and Brueckner— and Sham^ giving the result 
Z = 1.1978. The positive sign was unexpected, and gave 
a correction to the LDA that went in the wrong direc- 
tion, setting back the application of gradient corrections 
for years. 

The first step in the resolution of this puzzle came 
with Rasolt and Geldart's breakdown of Z according to 
the process involved, dividing Z = Z a (, + Z c , accord- 
ing to the processes depicted in Fig. la,b,c of Ref. [48] | 
and Fig . 3a,b,c of Ref. and also Fig. 4a,b,c of 

Ref. [46| where the same results were obtained by a dif- 
ferent method. The inset c in each of these figures shows 
what is known as the fluctuation diagram. It is the cause 
of the unphysical positive contribution to Z: one finds 
Z ab = -0.8491 and Z c = 2.0470. 

It would be another decade, however, before the im- 
plications of the fluctuation diagram would become fully 
apparent. This began with a seminal paper by Maggs and 
Ashcroftf^ who showed that the diagram was associated 
with the van der Waals interaction, followed by further 
development ) 47 ' 54 and culminating 55 with the use of the 
fluctuation diagram to obtain the London formula for 
the van der Waals interaction between two atoms. The 
reason for the failure of the gradient expansion was now 
apparent: the calculation of Z c represented an attempt 
to expand the long-range vdW interaction in powers of 
density gradients, an enterprise that was in retrospect 
obviously doomed to fail. 

Based on these facts, we write for use as the second 
term on the right side of Eq. (12) of Ref. [Ic| 



^grad — Z ai 



e 2 kp 



Vn 



127T \2kpn 



(B2) 



The van der Waals contribution is already included in 
another way, and it would be double counting to add the 
gradient expansion to it in addition. This is a matter 
of principle, but convenient indeed, because it avoids the 
use of Z c , which fails to be a valid approximation. This 
explains the use of Z a \, rather than Z in Eq. (JT)). 



In principle, Z a t> is a not a constant, but rather a func- 
tion of electronic density. There is no published data on 
this dependence. In Fig. 2 of Ref. [45| are shown two 
calculations that give the density dependence of the full 
Z, one of which was obtained earlie r 48 ' 49 [Z is equal to 
the ordinate of that figure times a numerical constant]. 
When the calculations of Ref. [45[ were made, the den- 
sity dependence of Z a b and Z c were calculated separately, 
but the significance of the individual quantities was not 
known at that time, and individual results were not pre- 
sented. However, the common author of this paper and 
Ref. [45| , who did this part of the calculation, attests that 
the principal density dependence shown in the Fig. 2 of 
this reference comes from that of Z c , and that the density 
dependence of Z a i, was weak. We thus feel comfortable 
using the high density value of Z a b = —0.8491. 



APPENDIX C: EVALUATION OF THE 
GRADIENT 

In the following, we will focus on the term V • [. . . ] in 
Eq. (Til)]) . It is straight forward to see that 

V- [^(r,r')i? rr 's(r)] = (V^(r, r'))R rT , • s(r)+ 
+0«,(r,r , )(ViJ„/)-8(r) + d (r,r / )iJ„/V-B(r) . 

The gradient of R rr i is the unit vector R rr <. For the 
gradient of </>d(r, r') we abbreviate higher derivatives of 
4>(d, d!) in the same fashion as above, 

d 2 ^(d,d>) d 2 ^d>) 

cj)dd{d, d ) and = (f> dd >{d, d ) . 



dddd 



dddd' 



Inserting the definition of d(r, r') and d'(r, r') from 
Eq. ©, the gradient of <j> d (r, r ') can then be written as 

V^(r,r') = <t> dd (r,r')Vd(r,r') + 4> dd ,(r,r')Vd'(r,r') 
= <f) d d{r, r')(R rr /g (r) + i? rr 'Vg (r)) 

+ fod' (r,r')Rrr'«o(r') . (CI) 

If we collect corresponding terms, the complete gradient 
in Eq. (fT5|) is then easily seen to be 



6 d (r,r')i2 Pr 's(r) = R rr ' ■ s(r) x 

</>d(r,r') + c/W(r,r')<7o(r)i?rr' + <?W( r > r')q (r')R r 
+ <Mr, r')^rr'S(r) • X7q (r) + &,(r, r')i? rr ' V • s(r) 



Expression for <p d (d, d'), 4>dd{d, d'), and <fidd'(d, d') are dis- 
cussed in Appendix lAl 
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